home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / unimin.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  9.8 KB  |  257 lines  |  [TEXT/????]
  1. ;;; modified by mh 12/6/87
  2. ;;; $Header: unimin.scm,v 1.4 88/01/26 07:45:38 GMT gjs Exp $
  3. ;;; 7/5/87 UNIMIN.SCM -- first compilation of routines for univariate
  4. ;;;             optimization.
  5.  
  6. (if-mit
  7.  (declare (usual-integrations = + - * /
  8.                  zero? 1+ -1+
  9.                  ;; truncate round floor ceiling
  10.                  sqrt exp log sin cos)))
  11.  
  12. ;;; Golden Section search, with supplied convergence test procedure,
  13. ;;; GOOD-ENUF?, which is a predicate accepting seven arguments: a, minx,
  14. ;;; b, fa, fminx, fb, count.
  15.  
  16. (define (golden-section-min f a b good-enuf?)
  17.   (define g1 (/ (- (sqrt 5) 1) 2))
  18.   (define g2 (- 1 g1))
  19.   (define (new-left lo-x hi-x) 
  20.     (+ (* g2 hi-x) (* g1 lo-x)))
  21.   (define (new-right lo-x hi-x)
  22.     (+ (* g2 lo-x) (* g1 hi-x)))
  23.   (define (best-of a b c fa fb fc)
  24.     (if (< fa (min fb fc))
  25.         (list a fa)
  26.         (if (<= fb (min fa fc))
  27.             (list b fb)
  28.             (list c fc))))
  29.   (define (lp x0 x1 x2 x3 f0 f1 f2 f3 count)
  30.     (if (< f1 f2)
  31.         (if (good-enuf? x0 x1 x2 f0 f1 f2 count)
  32.             (append (best-of x0 x1 x2 f0 f1 f2) (list count))
  33.         (let ((nx1 (new-left x0 x2)))
  34.           (lp x0 nx1 x1 x2 f0 (f nx1) f1 f2 (1+ count))))
  35.     (if (good-enuf? x1 x2 x3 f1 f2 f3 count)
  36.         (append (best-of x1 x2 x3 f1 f2 f3) (list count))
  37.         (let ((nx2 (new-right x1 x3)))
  38.           (lp x1 x2 nx2 x3 f1 f2 (f nx2) f3 (1+ count))))))
  39.   (let ((x1 (new-left a b)) (x2 (new-right a b)))
  40.     (let ((fa (f a)) (fx1 (f x1)) (fx2 (f x2)) (fb (f b)))
  41.       (lp a x1 x2 b fa fx1 fx2 fb 0))))
  42.  
  43.  
  44. ;;; For convenience, we also provide the sister-procedure for finding
  45. ;;; the maximum of a unimodal function.
  46.  
  47. (define (golden-section-max f a b good-enuf?)
  48.   (let* ((-f (lambda (x) (- (f x))))
  49.          (result (golden-section-min -f a b good-enuf?))
  50.          (x (car result)) (-fx (cadr result)) (count (caddr result)))
  51.     (list x (- -fx) count)))
  52.  
  53.  
  54. ;;; The following procedure allows keyword specification of typical
  55. ;;; convergence criteria.
  56.  
  57. (define (gsmin f a b . params)
  58.   (let ((ok? 
  59.            (if (null? params)
  60.                (lambda (a minx b fa fminx fb count)
  61.                  (close-enuf? (max fa fb) fminx (* 10 *machine-epsilon*)))
  62.                (let ((type (car params))
  63.                (value (cadr params)))
  64.                  (cond ((eq? type 'function-tol)
  65.                         (lambda (a minx b fa fminx fb count)
  66.                           (close-enuf? (max fa fb) fminx value)))
  67.                 ((eq? type 'arg-tol)
  68.                 (lambda (a minx b fa fminx fb count)
  69.                           (close-enuf? a b value)))
  70.                ((eq? type 'count)
  71.                 (lambda (a minx b fa fminx fb count)
  72.                   (>= count value))))))))
  73.     (golden-section-min f a b ok?)))
  74.  
  75. (define (gsmax f a b . params)
  76.   (let ((ok? 
  77.            (if (null? params)
  78.                (lambda (a minx b fa fminx fb count)
  79.                  (close-enuf? (max fa fb) fminx (* 10 *machine-epsilon*)))
  80.                (let ((type (car params))
  81.                (value (cadr params)))
  82.                  (cond ((eq? type 'function-tol)
  83.                         (lambda (a minx b fa fminx fb count)
  84.                           (close-enuf? (max fa fb) fminx value)))
  85.                 ((eq? type 'arg-tol)
  86.                 (lambda (a minx b fa fminx fb count)
  87.                           (close-enuf? a b value)))
  88.                ((eq? type 'count)
  89.                 (lambda (a minx b fa fminx fb count)
  90.                   (>= count value))))))))
  91.     (golden-section-max f a b ok?)))
  92.  
  93.  
  94. ;;; Brent's algorithm for univariate minimization -- transcribed from
  95. ;;; pages 79-80 of his book "Algorithms for Minimization Without Derivatives"
  96.  
  97. (define small-bugger-factor 1e-14)
  98.  
  99. (define (brent-min f a b eps)
  100.   (let* ((g (/ (- 3 (sqrt 5)) 2))
  101.          (a a)
  102.          (b b)
  103.          (x (+ a (* g (- b a))))
  104.          (fx (f x))
  105.          (w x) (fw fx) (v x) (fv fx)
  106.          (d 0) (e 0) (old-e 0) (p 0) (q 0)
  107.          (u 0) (fu 0))
  108.     (let loop ((count 0))
  109.       (let* ((tol (+ (* eps (abs x)) small-bugger-factor))
  110.              (2tol (* 2 tol))
  111.              (m (/ (+ a b) 2)))
  112.         ;; test for convergence
  113.         (if (< (max (- x a) (- b x)) 2tol)
  114.             (list x fx count)
  115.             (begin
  116.               (if (> (abs e) tol)
  117.                   (let* ((t1 (* (- x w) (- fx fv)))
  118.                          (t2 (* (- x v) (- fx fw)))
  119.                          (t3 (- (* (- x v) t2) (* (- x w) t1)))
  120.                          (t4 (* 2 (- t2 t1))))
  121.                     (set! p (if (positive? t4) (- t3) t3))
  122.                     (set! q (abs t4))
  123.                     (set! old-e e)
  124.                     (set! e d)))
  125.               (if (and (< (abs p) (* 0.5 q old-e))
  126.                        (> p (* q (- a x)))
  127.                        (< p (* q (- b x))))
  128.                   ;; parabolic step
  129.                   (begin (set! d (/ p q))
  130.                          (set! u (+ x d))
  131.                          (if (< (min (- u a) (- b u)) 2tol)
  132.                              (set! d (if (< x m) tol (- tol)))))
  133.                   ;;else, golden section step
  134.                   (begin (set! e (if (< x m) (- b x) (- a x)))
  135.                          (set! d (* g e))))
  136.               (set! u (+ x (if (> (abs d) tol) 
  137.                                d
  138.                                (if (positive? d) tol (- tol)))))
  139.               (set! fu (f u))
  140.               (if (<= fu fx)
  141.                   (begin (if (< u x) (set! b x) (set! a x))
  142.                          (set! v w) (set! fv fw)
  143.                          (set! w x) (set! fw fx)
  144.                          (set! x u) (set! fx fu))
  145.                   (begin (if (< u x) (set! a u) (set! b u))
  146.                          (if (or (<= fu fw) (= w x))
  147.                              (begin (set! v w) (set! fv fw)
  148.                                     (set! w u) (set! fw fu))
  149.                              (if (or (<= fu fv) (= v x) (= v w))
  150.                                  (begin (set! v u) (set! fv fu))))))
  151.               (loop (+ count 1))))))))
  152.  
  153. (define (brent-max f a b eps)
  154.   (define (-f x) (- (f x)))
  155.   (let ((result (brent-min -f a b eps)))
  156.     (list (car result) (- (cadr result)) (caddr result))))
  157.  
  158.                            
  159. ;;; Given a function f, a starting point and a step size, try to bracket
  160. ;;; a local extremum for f. Return a list (retcode a b c fa fb fc iter-count)
  161. ;;; where a < b < c, and fa, fb, fc are the function values at these
  162. ;;; points. In the case of a minimum, fb <= (min fa fc); the opposite
  163. ;;; inequality holds in the case of a maximum. iter-count is the number
  164. ;;; of function evaluations required. retcode is 'okay if the search
  165. ;;; succeeded, or 'maxcount if it was abandoned.
  166.  
  167. (define (bracket-min f start step max-tries)
  168.   (define (reorder a b c fa fb fc) ;return with a < c
  169.     (if (< a c)
  170.         (list a b c fa fb fc)
  171.         (list c b a fc fb fa)))
  172.   (define (test-it a b c fa fb fc count) ;assumes b between a and c
  173.     (if (> count max-tries)
  174.         `(maxcount ,@(reorder a b c fa fb fc) ,max-tries)
  175.         (if (<= fb (min fa fc))
  176.             `(okay ,@(reorder a b c fa fb fc) ,count)
  177.             (let ((d (+ c (- c a))))
  178.               (test-it b c d fb fc (f d) (+ count 1))))))
  179.   (let ((a start) (b (+ start step)))
  180.     (let ((fa (f a)) (fb (f b)))
  181.       (if (> fa fb)
  182.           (let ((c (+ b (- b a))))
  183.             (test-it a b c fa fb (f c) 0))
  184.           (let ((c (+ a (- a b))))
  185.             (test-it b a c fb fa (f c) 0))))))
  186.  
  187. (define (bracket-max f start step . max-tries)
  188.   (define (-f x) (- (f x)))
  189.   (apply bracket-min (append (list -f start step) max-tries)))
  190.  
  191.  
  192. ;;; Given a function f on [a, b] and N > 0, examine f at the endpoints
  193. ;;; a, b, and at N equally-separated interior points. From this form a
  194. ;;; list of brackets (p q) in each of which a local maximum is trapped. 
  195. ;;; Then apply Golden Section to all these brackets and return a list of
  196. ;;; pairs (x fx) representing the local maxima.
  197.  
  198. (define (local-maxima f a b n ftol)
  199.   (let* ((h (/ (- b a) (+ n 1)))
  200.          (xlist (generate-list 
  201.                   (+ n 2) 
  202.                   (lambda (i) (if (= i (+ n 1)) b (+ a (* i h))))))
  203.          (flist (map f xlist))
  204.          (xi (lambda(i) (list-ref xlist i)))
  205.          (fi (lambda(i) (list-ref flist i)))
  206.          (brack1 (if (> (fi 0) (fi 1))
  207.                      (list (list (xi 0) (xi 1)))
  208.                      '()))
  209.          (brack2 (if (> (fi (+ n 1)) (fi n))
  210.                      (cons (list (xi n) (xi (+ n 1))) brack1)
  211.                      brack1))
  212.          (bracketlist 
  213.            (let loop ((i 1) (b brack2))
  214.              (if (> i n)
  215.                  b
  216.                  (if (and (<= (fi (- i 1)) (fi i)) 
  217.                           (>= (fi i) (fi (+ i 1))))
  218.                      (loop (+ i 1) (cons (list (xi (- i 1)) 
  219.                                                (xi (+ i 1))) b))
  220.                      (loop (+ i 1) b)))))
  221.          (locmax (lambda (int) (gsmax f (car int) (cadr int) 
  222.                                       'function-tol ftol))))
  223.     (map locmax bracketlist)))
  224.                    
  225. (define (local-minima f a b n ftol)
  226.   (let* ((g (lambda (x) (- (f x))))
  227.          (result (local-maxima g a b n ftol))
  228.          (flip (lambda (r) (list (car r) (- (cadr r)) (caddr r)))))
  229.     (map flip result)))
  230.        
  231.  
  232. (define (estimate-global-max f a b n ftol)
  233.   (let ((local-maxs (local-maxima f a b n ftol)))
  234.     (let loop ((best-so-far (car local-maxs)) 
  235.                (unexamined (cdr local-maxs)))
  236.       (if (null? unexamined)
  237.           best-so-far
  238.           (let ((next (car unexamined)))
  239.             (if (> (cadr next) (cadr best-so-far))
  240.                 (loop next (cdr unexamined))
  241.                 (loop best-so-far (cdr unexamined))))))))
  242.  
  243. (define (estimate-global-min f a b n ftol)
  244.   (let ((local-mins (local-minima f a b n ftol)))
  245.     (let loop ((best-so-far (car local-mins)) 
  246.                (unexamined (cdr local-mins)))
  247.       (if (null? unexamined)
  248.           best-so-far
  249.           (let ((next (car unexamined)))
  250.             (if (< (cadr next) (cadr best-so-far))
  251.                 (loop next (cdr unexamined))
  252.                 (loop best-so-far (cdr unexamined))))))))
  253.  
  254.  
  255.     
  256.  
  257.